N6831516 


FORTRAN SUBPROGRAMS FOR COMPLETE ELLIPTIC INTEGRALS 


by 


F. I. Zonis 

RCA laboratories 
Princeton., New Jersey 


The Research discussed in this paper was partially sponsored 
by The National Aeronautics and Space Administration, Ames 
Research Center, California, under Contract Number NAS2-3772. 





GENERAL DISCLAIMER 


This document may have problems that one or more of the following disclaimer 

statements refer to: 


This document has been reproduced from the best copy furnished by the 
sponsoring agency. It is being released in the interest of making 
available as much information as possible. 

This document may contain data which exceeds the sheet parameters. It 
was furnished in this condition by the sponsoring agency and is the best 
copy available. 

This document may contain tone-on-tone or color graphs, charts and/or 
pictures which have been reproduced in black and white. 

The document is paginated as submitted by the original source. 

Portions of this document are not frilly legible due to the historical nature 
of some of the material. However, it is the best reproduction available 
from the original submission. 



RCA LABORATORIES 


Radio Corporation of America 
Princeton, Kev Jersey 

FORTRAN SUBPROGRAMS FOR COMPLETE ELLIPTIC INTEGRALS 

by 
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ABSTRACT 

Fortran IT subprograms have been developed for evaluating the complete 
elliptic i ntegrals of the first and second kinds on the RCA 601 computer. 
These subprograms exhibit an error no greater than 2 (10) 
entire range of definition. 


over the 




DEFINITIONS 


The complete elliptic integral of the first kind is defined as 


K(m) = J [(l - t 2 ) (l - mt 2 )] 
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and the complete elliptic integral of the second kind is defined as 
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In the above expressions m is the parameter of the integrals. This quan- 
tity is related to the modulus k and the modular angle Ct by the relations 


i 2 . 2 

m = k - sm a. 
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One may also define the complementary parameter and compl cmentar y 
modulus k* through the relations 


m^ - 1 - m 
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= (k* ) = cos a 
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The user is cautioned that many texts define elliptic integrals in terms 
of the modulus k. 


PROGRAMMING METHOD 


The complete elliptic, integrals are evaluated using the polynomial approx 
imations given by Equations 17.3.34 and 17.3.36 of the Handbook of Mathe - 
matical Functions : ^ ^ 


17.3.34 
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Listings of the subprograms are given in Figures 1 and 2. In these sub- 
programs the Fortran floating point variable E plays the role of the com- 
plementary parameter The natural logarithm subroutine LOG(B) in the 

RCA 601 Fortran II package normally causes a loss of significant figure 

-2 

accuracy halt if 1-B <10 . In the present case^ however } the log- 

arithm will maintain sufficient accuracy when used in an expression of 
the form 


£ a . nr - X n (nr ) £ b. im 

x 1 v V . i 1 

i r 


with 


£ a . m 1 & £ b . m 1 ‘ 
l 1 l 1 

so that the result will be accurate to full significance. Thus the loss 
of accuracy halt was inhibited in these subprograms by using the special 
call BYPASS (LOG jB) . 


CALLING PROCEDURE 

The complete elliptic integrals of the first and second kind are evaluated 
by placing the terms ELK(B) and ELE (B ) y respectively ; in any floating- 
point Fortran arithmetic expression. Note that the calls are in terms of 
the complementary parameter B = m^. This was done^ on the suggestion of 
R. W. Klopf enstein^ to avoid the loss of accuracy in the machine computa- 
tion of 

B = 1 - (1 - B) (5) 

when B is known to full accuracy. (Note in Figures 1 and 2 that the in- 
tegrals are evaluated in terms of B.) 

EXAMPLE: To evaluate the expression 

C = K (ra = A) -I- E (n^ = B) + 2 K (a = D°) (6a) 

1 

with A, B; and D known, one could write the Fortran statement 
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C = ELK (1.0 - A) + ELE (B) + 2.0 * ELK (COS DF (D) ** 2) 


(6b) 


ERROR STOPS AND SPECIAL CONDITION 


B must satisfy 0 < B < 1 for ELK (B) and 0 < B < 1 for ELE (B) . For B 
outside these limits an error message is printed giving the value of B. 
The job is then terminated with a dump. 

For B = 0 in ELE (B) the polynomial approximation is bypassed, ELE is set 
equal to 1.0, and control is returned to the calling program. 

ACCURACY 

Provided B is given to full (i.e., 9 significant figures) accuracy, ELE 

- 8 

and ELK will exhibit an error of no more than 2 (10) > that is, 2 units 

in the ninth significant figure* 

In testing the ELK program, two situations were encountered where poorer 
accuracy was obtained. First, as would be expected, calls of the form 
ELK (l.O - (1.0 - B)j for small B yielded results which were accurate to 
the same number of significant figures as ^1.0 - (1.0 - B)^ . Thus, the 
above form of the call should be avoided when the complementary parameter 
B is available. Second, calls of the form ELK (COS DF (D) ** 2) lost 
several significant figures of accuracy when D was very close to 90 . 

This is due to a loss of significance, in COS DF at these values, as is 
illustrated below. 

ELE is not subject to the same loss of significance since in this program 
there is no constant term in polynomial which multiplies the logarithm. 
[See 17.3.36 of Reference (1).] 


- 
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TESTING 


Bo tli ELK and ELE were evaluated for the following values of m and a 
m = 10 5 , 10 \ 10 3 , 10 2 , 0.1(0.05)0,9(0.01)0.99, 0.999, 0.9999, 0.99999 
a = 0 (10) SO (1) 89 (0.1) 89.9 degrees 
The results are shown in Figures 3 through 8* 

The computed values we re compared with values taken from Tables 17.1 and 
17.2 of Reference (1) or values computed by a special double precision 
program described below. . The A columns following the computed values of 
ELK = FIRST and ELE = SECOND in Figures 3 through 8 give 

& - (e - E \ (10) 8 (7) 

\ comp. exacty v ' 

where E is the value of the elliptic integral. 

It can be seen that both subprograms maintain the specified accuracy over 
the entire range of m. However, the elliptic integral of the first kind 
looses some significance for a close to 90°. That this is due to a loss 
of significance in the cosine evaluation can be seen in Figures 6 through 
4 8 where the column headed gives the error in the ninth significant 

figure of the computed value of the cosine. Note in particular the loss 
of significant figures in Figure 8. 

Figure 9 shows the results of the evaluation of ELE (O' = 90°) and also 
shows the error message printed out when ELK (90°) was called. 


SPECIAL TEST PROGRAM 

In order to obtain accurate values of the elliptic integrals outside the 
range covered in the tables, a special test program was written in double 
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precision for the 70/45 Phase I Basic Time Sharing System. This program 
is shown as Figure 10. 

This program uses Equation (8) to obtain three stages of reduction of the 
parameter m: 


where 


and 


m. , 
l + l 


/ 1 - /l~- nn 


1 1 + f\. + in. 
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m. 
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( 8 ) 


( 9 ) 

( 10 ) 


Next, K(m ) and E(m ) are evaluated using the series expansions 773.2 and 
8 ( 2 ) 

774.2 of Dwight. ' (The reduction on m assures that these expansions 
are rapidly convergent, even for m very close to 1.) Finally, Equations 
17.3.29 and 17.3.30 of Reference (1) are applied three times to obtain 
K(m) and E(m), respectively. 

Figure 11 shows the results obtained from this program for selected val- 
ues of m. Comparison with tabulated values show the results are accurate 
to 13 significant figures. 

To obtain Figure 12, the statement 


35 M = SIN(1.5707963267948966*M/90)**2 

was inserted between statements 30 and 40 in the test program. Thus, in 
this table M represents the modular angle (in degrees). Better than 9 
significant figure accuracy was obtained over the entire range. 


_ 6 _ 



The coding procedure used in the test program was not used in the subpro- 
grams since it requires more code and takes longer to execute. 
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